COMBINED PRECONDITIONING WITH APPLICATIONS IN 
RESERVOIR SIMULATION 
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Abstract. Wc develop a simple algorithmic framework to solve large-scale symmetric positive 
definite linear systems. At its core, the framework relies on two components: (1) a norm-convergent 
iterative method (i.e. smoother) and (2) a preconditioner. The resulting preconditioner, which we 
refer to as a combined preconditioner, is much more robust and efficient than the iterative method 
and preconditioner when used in Krylov subspace methods. We prove that the combined precondi- 
tioner is positive definite and show estimates on the condition number of the preconditioned system. 
We combine an algebraic multigrid method and an incomplete factorization preconditioner to test 
the proposed framework on problems in petroleum reservoir simulation. Our numerical experiments 
demonstrate noticeable speed-up when we compare our combined method with the standalone alge- 
braic multigrid method or the incomplete factorization preconditioner. 
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1. Introduction. We consider a sparse linear system of equations that arises 
from the discretizations of the elliptic partial differential equations (PDEs) used to 
simulate the physical processes in petroleum reservoirs. The Petroleum Reservoir 
Simulation (PRS) is a tool for predicting hydrocarbon reservoir performance under 
various operating regimes. PRS helps engineers to obtain information pertaining to 
the processes that take place within oil reservoirs - information that can be used to 
maximize recovery and minimize environmental damage. The crucial aspect of PRS is 
its ability to solve large-scale discretized PDEs, which are strongly coupled, indefinite 
and often non-symmetric. Solving these linear systems consumes most of the compu- 
tational time in all modern reservoir simulators - more than 75% of the computational 
time in general. Furthermore, the demand for more accurate simulations has led to 
larger discrete reservoir models. And, this increase in model size results in larger 
linear systems, that are more difficult, or even impossible to solve in an acceptable 
amount of time using standard direct or standard iterative solvers. 

Over the last 30 years, incomplete LU factorization (ILU) has become one of 
the most commonly used methods for solving large sparse linear equations arising in 
PRS. First developed in the 1960s [HJ H7]> ILU methods provide an approximation 
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of the exact LU factorization (computed via Gaussian elimination) by specifying the 
sparsity of the L and U factors. These methods need not be convergent when used 
in linear iterative procedures, but they do provide preconditioners that are used in 
Krylov subspace methods. Thought to have introduced the word preconditioning, 
Evans [11] may also have been the first to use ILU as a preconditioner. As noted, due 
to their simplicity, ILU methods are of particular interest to researchers in the field of 
reservoir simulation. However, when ILU preconditioners are applied to problems in 
PRS, their performance usually deteriorates as the number of grid-blocks increases. 

To solve discretized scalar elliptic PDEs (Poisson-like PDEs), multigrid (MG) 
methods are efficient and provide a solution in optimal time and memory complexity 
in cases that allow standard coarse spaces to be used (see [131 132 [21 [25] and references 
therein for details). However, the multigrid methods with standard coarse spaces make 
extensive use of the analytic information from the discretized equation and geometric 
information explicitly related to the discretization. This makes such methods difficult 
to use, such that in practice more sophisticated methods, such as algebraic multigrid 
(AMG) methods are preferred. There are many different types of AMG methods - the 
classical AMG ([21H]), smoothed aggregation AMG (HE GO), AMGe (QHCEB]), etc. 
However, despite their differences, they generally do not require geometric informa- 
tion from the grids. Interested readers can refer to [201 El El EH EH] and references 
therein for details on AMG methods. Due to their efficiency, scalability, and appli- 
cability, AMG methods have become increasingly popular in practice (see e.g., [24] ) 
including in modern petroleum reservoir simulations [23 , 3 EH EU [33] ■ However, the 
performance and efficiency of MG methods may degenerate as the physical and geo- 
metric properties of the problems become more complex. In such circumstances, there 
is an increase in AMG-setup time (the time needed to construct coarse spaces and 
operators) and superfluous fill-in in the coarse grid operators, which makes applying 
relaxation more expensive on coarser levels. 

In this paper, we provide a simple and transparent framework for combining pre- 
conditioners like ILU or Jacobi or additive Schwarz preconditioners with AMG or 
another norm-convergent iterative method for large-scale symmetric positive definite 
linear system of equaitons. In the combined preconditioners we propose, the com- 
ponent provided by the norm-convergent iterative method need not be very effective 
when used alone. However, as we shall see, when it is paired with an existing pre- 
conditioner (such as an ILU preconditioner), the result is an efficient method. We 
show that the combined preconditioner is positive definite (SPD) under some mild as- 
sumptions. This guarantees the convergence of the conjugate gradient (CG) method 
for the preconditioned system. We apply the combined preconditioner to petroleum 
reservoir simulations that involve highly heterogeneous media and unstructured grids 
in three spatial dimensions (with faults and pinch-outs). The numerical results show 
an improved performance that justifies using the combined preconditioners in PRS. 

The rest of the paper is organized as follows. The algorithmic framework is 
introduced and analyzed in Section 2. As an example in Section 3, we formally 
present the combined preconditioner by describing its essential components, e.g., an 
ILU preconditioner and an AMG method. In Section 4, we recall the mixed-hybrid 
finite clement method used to provide a discretization of an important component 
of reservoir simulation - an elliptic PDE with heterogenous coefficients. Numerical 
results showing the improvement in performance are presented in Section 5, and some 
concluding remarks are given in Section 6. 
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2. Definition of the combined preconditioner. Consider the linear system 
(2.1) Au = f, 

with an SPD coefficient matrix A. Let (•, •) be an inner product on a finite-dimensional 
Hilbert space V; its induced norm is || • ||. The adjoint of A with respect to (•,•), 
denoted by A T , is defined by (Au, v) = (u, A T v) for all u,v eV. A is, SPD if A T = A 
and (Av,v) > for all v <E V^\{0}. As A is SPD with respect to (•, •), the bilinear 
form (A-, ■) defines an inner product on V, denoted by (•, -)a, and the induced norm 
of A is denoted by | • || a- 

In what follows, we also need an operator S, referred to as a smoother. In some 
of the results, we distinguish two cases: (1) norm-convergent smoother, i.e., (I — SA) 
is a contraction in || • H^-norm; and (2) more generally, a non-expansive smoother, 
i.e., (I — SA) is a non-expansive operator in || ■ \\a norm. Note that norm-convergent 
implies non-expansive. 

Next, in Algorithm |2.1| we combine a non-expansive iterative method S with an 
existing SPD preconditioner B, and as a result we get a preconditioner B co . 

Algorithm 2.1 

Given S, B, and u k '° = u k , the new iterate u k+1 is obtained by the following steps: 

(1) u^ 1 = u k >° + S{f - Au k -°); (2) u k ' 2 = it**' 1 + B(f - Au k ^); 
(3) u k+1 = u k - 2 + S T (f - Au k < 2 ). 



It is easy to see that Algorithm |2.1| defines a linear operator B co and that 
(2.2) / - B co A = (I- S T A)(I - BA){I - SA). 



From \2.2\ , it follows that 

B co = S + (B - S T AB - BAS + S T ABAS) 
(2.3) =S+(I— S T A)B(I — AS). 

Here S is the symmetrization of S, defined as 

I -SA = {I - S T A)(I - SA), S = S + S T - S T AS. 

We now explain the properties pertaining to the notions of norm-convergent and 
non-expansive smoothers, which we refer to in the subsequent sections. 

• Norm-convergent smoother. In this case, (/ — SA) is a contraction in 
|| ■ Hyi-norm; i.e., there exists p e [0, 1), such that for all v e V we have 

(2-4) \\(I-SA)v\\l<p\\v\\ 2 A . 

Note that under this assumption, S is invertible. Indeed, we have, 

\\sav\\ a = (/- sa))v\\ a > NU - IK' - SA)v\\ A > (1 - Vp)IMU, 

which means SA is coercive and implies that S is invertible. Recall also 
that (2.4) holds if and only if the symmetrization S is SPD; i.e., 

\\I-SA\\ A <p, if and only if (Sv,v) > 0, for all veV\{0}. 



4 



• Non-expansive smoother. This is a more general case in which we assume 
that 

(2.5) \\(I-SA)v\\ A <\\v\\ A VveV. 

One important difference between these two cases is this: with the non- 
expansive case, S can be positive semi-definite, whereas in the norm-convergent 
case, S is positive definite. 
The following result shows that the operator J5 co is SPD and, therefore, can be em- 
ployed as a preconditioner for the preconditioned conjugate gradient (PCG) method. 
A preliminary version of the following theorem can be found in [33]. 

Theorem 2.1. Assume that S : V — > V is a non- expansive smoother. Moreover, 



the operator B : V —> V is SPD. Then, B co defined in (2.2) is SPD 



Proof. For any x € V, x ^ 0, we define x := (I — SA)x. Thus, we obtain 

((/ - B co A)x, x) A = ((I - BA)(I - SA)x, (I - SA)x) A 

= ((I - BA)x, x) A = (x, x) A - (Ax, Ax) B 

where (•, -) B is an inner product on V defined as (•, -)b ■= (B-, •) since B is SPD. 

Further, if x = 0, then ((/ — B co A)x, x) A = < ||x||^. On the other hand, for 
x 7^ 0, we obtain from (Ax, Ax) b > that 

\\x\\a - (B co Ax,x) A = ((I — B co A)x, x)a < \\x\\ A < \\x\\ A , (B co Ax,Ax) > 0. 

As x € V is an arbitrary nonzero element, and A is non-singular, we conclude that 



B co is SPD. Finally, in view of (2.3) and because both A and B are symmetric, it is 
easy to see that B co is symmetric. This completes the proof. □ 

Remark 2.1. Note that the assumption in the theorem is that the error transfer 
for the convergent iteration I — SA is non-expansive in the A norm and does not 
require I — SA to be a contraction. 
Remark 2.2. Theorem : 



2.1 



shows that the operator B co is SPD. Hence, we can use 
the PCG method to solve (2.1 ) with B co as a preconditioner, and in exact arithmetic, 
PCG method will always be convergent. Algorithm 2.1 is relatively simple, but the 
order in which S and B are applied is crucial to the positive definiteness property of 
B co . For example, consider the operator B defined by 

(2.6) I-BA=(I- BA) (I - SA)(I - BA) . 

In this case, it is not true in general that B is positive definite. For example, if we 
use the ILU method to define a preconditioner B, finding the right scaling to ensure 
the positive definiteness of B could be a difficult task, and with indefinite B the PCG 
will not converge. 



Remark 2.3. Note that directly using Algorithm 2.1 as an iterative method may 
not result in a convergent method. This is because though we assume that B is only 
SPD, we do not assume that I — BA is a contraction. The contraction property of 
/ — BA (and also of J — B co A) is that any eigenvalue of BA satisfies < X(BA) < 
uj < 2. However, we do not assume that such a contraction property holds for B. 

2.1. Other derivations of the combined preconditioner. Here, we present 
two distinct points of view that can lead to algorithms such as Algorithm |2.1| One 
is from the point of view of the fictitious or auxiliary space techniques, developed 
in [TH] and [32], and another is from the point of view of the block factorization 
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techniques's, developed in (28j Chapter 5]. Below, we derive Algorithm 2.1 using 
these techniques. Of course, all three derivations lead to one and the same method 
and they are equivalent to each other. However, they present different points of view 
and help us understand different aspects of the method. 

2.1.1. Derivation via the auxiliary space method. An additive version of 
the Algorithm |2.1| was derived via fictitious or auxiliary space techniques developed 



in [THj and [55]. Equivalently, Algorithm 2.1 can be viewed as a successive or multi 



plicative version of the auxiliary space additive preconditioner described below. Here, 
we use only one auxiliary space, which will turn out to be the same as V, but with a 
different inner product. We define 

(2.7) V = VxWi, 

where W± is an auxiliary (Hilbert) space with an inner product 

ai{-, •) = (-, -) Al - 

We also introduce an operator Hi : W\ h-> V in order to define the additive precondi- 
tioner, and define the latter as follows: 

(2.8) b = s + niA^nf . 



A distinctive feature of the auxiliary space method is the presence of V in (2.7 1 
as a component of V and the presence of the symmetric positive definite operator 
S : V i— !► V. S is assumed to be SPD in order to guarantee that B is SPD and can be 
applied as a preconditioner for the PCG method. It can be proved that the condition 
number of the preconditioned system can be bounded as follows (see |32|): 

(2.9) k(BA) <<*(£+<*), 

if 

\\U lWl \\ 2 A < cjWw^, w 1 EW 1 , 

and 

(SAv,v) A < c 2 s (v,v) A , VveV. 

Moreover, for each v G V, there exists w\ € W\ such that 

v = Vq+H\Wi, v Q =v-TliWi, and (SAv ,v ) A + (w 1 ,w 1 ) Al <cl(v,v) A . 

Taking W\ = V but with the inner product defined by B, A^ 1 = B, and n = /, gives 
the additive version B of the preconditioner B co . Note that when W\ = V, B is SPD 
and S is the symmetrization of the smoother S, the preconditioner B is SPD for both 
the norm-convergent smoother and the non-expansive smoother. 

2.1.2. Derivation via block factorization. In this section, we present a deriva- 
tion following the lines in [28l Chapter 5] . Let us introduce the operator B : V x V i— > 
V X V in block factored form: 

;i -S T A\ (S 0\ ( I 

D = 



/ B \-AS I 



() 



Wc then set 



(2.10) B = (I, I)B ( 1 ) , B :V i->V. 



A straightforward calculation shows that B co and B are the same. Regarding the 
positive definiteness of B, note that for the norm-convergent smoother, it is immedi- 
ately evident that B is SPD, because both B and S are SPD. For the non-expansive 
smoother, the positive-definiteness of B does not immediately follow from the form 



of the preconditioner given in (2.10). The arguments of Theorem 2.1 (or similar) are 



needed to conclude that B is SPD. 

2.2. Effectiveness of B co . In this section, we show that the combined precon- 
ditioner, under suitable scaling assumptions performs no worse than its components. 
As the numerical tests show (see Section [5]), the combination of ILU (for B) and AMG 
(for S) performs significantly better than its components. Let us set 

(2.11) mi = X max (BA) and m = X min (BA). 

Without loss of generality (with proper scaling), we can assume that the precondi- 
tioner B is such that the following inequalities hold: 

(2.12) mi > 1 > m > 0. 

We now prove a result that compares k(B co A) with k(SA) and n(BA) under the 
assumptions that B and S are such that both ( |2.12 1 and ( |2.4[ ) are satisfied 



Theorem 2.2. If S is a norm-convergent smoother and B co is defined as in ( 2.2 ) 
then 

(2-13) K (B co A) < ^f-^ 1 , 

(1 - m )(l - p) + m 



(2.14) k(B co A) < k(BA), 



Furthermore, if S is such that (2.4 1 holds with p > 1 — m ™^, 1 > then 



(2.15) k(B co A) < k(SA). 



Proof. From the assumption stated in (2.12), we immediately conclude that B is 
SPD and k(BA) = mi /mo. By the definition of S, we have 

< ((/ - SA)w, w) A = ((I - SA)w, (I - SA)w) A = - SA)w\\\ < p\\w\\% 

where we have used the assumption of the convergence of S in the last inequality. By 
choosing v — Aw, we can obtain 



(2.16) 



(1 - p)(A~ 1 v, v) < (Sv, v) < (A~ 1 v, v). 
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On the other hand, as mi = \ maK (BA) and mo = A m i n (BA), we have 

fflo(A _1 !),n) < (Bv,v) < mi(A^ 1 v,v). 

By the definition of B co ( |2.3[ ), we have 

{B co v, v) = {Sv, v) + (B{I - AS)v, {I - AS)v) 

< (Sv,v) + m 1 (A- 1 (I - AS)v,(I - AS)v) 
= (1 — mi)(Sv, v) + mi(A^ 1 v, v). 

Similarly, we can derive that 

(B co v, v) = (Sv, v) + (B(I - AS)v, {I - AS)v) 

> {Sv,v) + m (A- 1 (I- AS)v,(I- AS)v) 
= (1 — m )(Sv, v) + m (A _1 w, v). 

As mi > 1 > mo, we have 

[(1 - m )(l - p) + m Q ](A- 1 v, v) < (B co v, v) < [(1 - m x )(l - p) + m^A^v, v); 
then the condition number of B co A is bounded by 

(1 - mi)(l - p) + mi 



If mi > 1, then 



and if 1 > mo > 0, then 



Note that 



n{B co A) < 



<B co A)< 

(1 - m )(l - p) + m 

(1 - mi)(l — p) + mi < mi, 

(1 - m )(l - p) + m > m . 
(1 — mi)(l — p) + mi mi 



< 



(1 - m )(l - p) +m m 



k(BA). 



Hence, the inequality (2.14| holds if mi > 1 > mo > 0. 
On the other hand, (2.151 follows from 

(2.17) 



(1 - mi)(l - p) +mi 1 



(1 - m )(l - p) + m l-p 



k(SA), 



where the last equality comes from (2.161. Note that 



(1 - mi)(l — p) + m 1 _ 1 + p(mi — 1) 
(1 - m )(l - p) + m (l-^ + pmo' 



mi —1 



Note that f±£ < f if £ < s. and a, b, c, d > 0. Therefore, (2. 17 1 holds if . . 



< 



i.e., /O > 1 



-. □ 



Remark 24. If =^ < j^, then ( |2.15[ ) becomes k(B co A) < k{SA). 
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Remark 2.5. If either S or B works well as a preconditioner itself, there is no 
need to use the more complicated B co . However, we are interested in cases in which 
neither S nor B works effectively alone. In these cases, we use B co to combine these 
two, and we use Theorem |2.2| to guarantees that B co will be a better preconditioner 
than either B or S under reasonable conditions. The condition 

p > 1 - 7710/(771! - 1) 

means that both S and B yield slow convergence as 

p > 1 - mo /(mi -l)wl - k(BA)^ 1 if m x » 1. 

Therefore, if B is not a good preconditioner, i.e., k(BA) is large, then p ss 1, which 
implies^ that S does not work well either. The new preconditioner B is no worse than 
B or S alone. In fact, the new preconditioner may be capable of performing much 
better than either B or S based on our numerical experiments (see Section 5b. 



3. A simple example: ILU+AMG. Algorithm |2. 1| provides an approach to 
combining ILU and AMG methods. In the combined use of ILU and AMG for the 



linear system (2.1 1, AMG serves as S and ILU serves as B. Next, we specify our 
choice of ILU and AMG for problems specific to reservoir simulation. However, we 
would like to emphasize that the combined preconditioner works for a wide range 
of iterative methods and preconditioners as long as they satisfy the assumptions in 
Theorem O 



3.1. Incomplete LU factorization. Incomplete LU factorizations compute a 
sparse lower triangular matrix L and a sparse upper triangular matrix U so that the 
residual matrix R = A — LU satisfies certain conditions. A general algorithm of ILU 
can be obtained by Gaussian elimination and dropping some of the elements in the 
off-diagonal positions. There are many variants of ILU preconditioners. They differ 
in terms of the rules that govern the prescribed fill-in of the factors during the ILU 
factorization procedure. For example, discarding the fill-in based on position gives the 
ILU(fc) method; discarding the fill-in based on the values of the corresponding entries 
in the L or U factors gives the threshold ILU method. There are also ILU methods 
for which the fill-in is managed based on a combination of positions and values or 
based on other dropping strategies. For details about the different ILU methods, we 
refer to monographs [IJE] and Saad [21]. Here we use the notation and terminology 
from Saad EU. 



We consider ILU(fc) (Algorithm 3.1), which to our knowledge was originally in- 
troduced for reservoir simulations in [30] . The ILU is based on the level of fill to 
determine the off-diagonal positions in the L and U factors where the entries fill-in 
will not be introduced (or as is often said, off-diagonal positions for which the fill-in 
entries are dropped). Next, we define level of fill, and the detailed algorithm is given 
in Algorithm |3.1| 

Definition 3.1 (Level of fill). The initial level of fill of the elements of a sparse 
matrix A defined as 

Cij = if a.ij ^ or i — j, otherwise Cij = oo. 

When an entry is updated in the factorization procedure ay := a,j — a^. * a^j, its 
level of fill is also updated by 



dj = min{Cij, C ik + C kj +!}■ 



9 



Remark 3.1. Although the level of fill depends on the location of the element, the 
rationale is that the level of fill should indicate the magnitude of the element. The 
higher the level of fill, the smaller the element. For example, when the level of fill of 
a,ij is k this means that the size is |aj,-| = 0(e k ) for some e < 1. 



Algorithm 3.1 ILU(fc) Method 

1. For all nonzero elements Ojj, define Cij = 

2. For i = 1, ■ ■ ■ n, Do 

For m = 1, • • • , i — 1 and C; Lm < k, Do 
Compute a im := a im /a mm 

Compute cii* := — ai m a mie , where a.^ is the i-th row of A 
Update the levels of fill of non-zero by Cij = minl/^-, Cik + Ckj + 1} 
End 

Replace any element in row i with Cij > k by zero 
End 



3.2. Algebraic multigrid methods. As an algebraic variant of MG methods, 
AMG methods are widely applicable and the focus of current intensive development. 
AMG methods have mesh-independent convergence rates and optimal computational 
complexity for a wide range of problems. 

Any AMG method consists of two phases: the SETUP phase and the SOLVE 
phase. In the SETUP phase, the intergrid operator Pi is constructed, and the coarse 
grid matrix is defined as 

Ai = PfA l+1 P u Z = L-1,...,0, 

and Al = A. In the SOLVE phase, the smoother Si and coarse-grid correction are 
applied recursively as shown in the following general V-cycle MG Algorithm |3.2| 



Algorithm 3.2 V-cycle for solving Aiui — ft, with an initial guess m° 

1. Pre-smoothing: u] = + Si(fi — Aiu®) 

2. Coarse-grid correction: 

(a) f l - l =Pl l {f l -A l u\) 



(b) If I = 1, eo = A 1 /o; else, apply Algorithm 3.2 for A;_ie/_i = with 
zero initial guess 

(c) uf = u\ + Pi-id-i 
Post-smoothing: uf = uf + Sf(fi ~ Aiuf) 



In AMG methods, Pi must be constructed based on algebraic principles, which 
presents certain challenges. The key to fast convergence is the complementary nature 
of operators Si and Pi. That is, errors not reduced by Si must be interpolated well 
by Pi. We choose the classical AMG as the iterative method S in our Algorithm |2.1| 
This method constructs the intergrid operator Pi in two steps. First, the classical 
C-F splitting method is used, and then the operator Pi is constructed by classical 
Ruge-Stiiben interpolation; see [2D] for details. 

3.3. Application of Algorithm |2.1|. In this paper, we consider SPD coeffi- 
cient matrices A only. In this case, ILU(fc) is replaced by an Incomplete Cholesky 
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(IC) factorization . It is easy to see that IC(fc) is SPD, which satisfies the assumption 
on B in Theorem 



2.1 



For AMG methods, we have A t = PfA l+1 P h I = L - 1, . . . , 0. 
A is SPD and Pi is constructed in the classical AMG method; therefore, it is guar- 
anteed that Ai is also SPD for I = 0,1, ... ,L. Hence, it is easy to see that using 
the standard Gauss-Seidel (GS) smoother on each level gives a convergent iterative 
method. According to Algorithm |3.2[ Bi, which stands for V-cycle MG on level I 
(I = 1, 2, . . . , L), can be defined recursively as 

/ - = (I- SfAi)(I - P l _ x B l _ 1 Pl_ 1 A l ){I - SiAi), 

with B = Aq 1 . As 

||7 - S l A l \\ Al < 1 and \\I - P t _ 1 B l _ 1 P?_ 1 A l \\ Al < 1, 

by induction assumption, we can show that 

iii-s^m, < i. 

Thus the classical AMG is convergent for SPD problems by mathematical induction. 



Hence, ILU(fc) and the classical AMG can be combined using Algorithm 2.1 to yield 
the following corollary: 

Corollary 3.2 (Symmetric Positive Definiteness of ILU-AMG). If we choose 



S from the classical AMG methods and B from ILU(k) in Algorithm 2.1 then the 
operator B co defined in (12. 21) is SPD. 



4. A model problem in reservoir simulation. Petroleum reservoir simula- 
tion provides information about processes that take place within oil reservoirs. It is, 
therefore, of great assistance in efforts to achieve optimal recovery. Modern reservoir 
simulation faces increasingly complex physical models and uses highly unstructured 
grids. This results in Jacobian systems that are more difficult at each step. In this 
section, we describe a model problem in PRS and a discretization method that has 
great promise for efficient resource recovery, which is used for numerical tests in the 
next section. 

Let f2 be a bounded domain in K 3 in our consideration of the following second- 
order elliptic problem: 

(4.1) - V • (aVp) + cp = / inO. 



In reservoir simulation, the model problem (4.1) usually occurs after temporal semi- 
discretization of the mathematical models that describe the multiphase flow in porous 
media. Such scalar problems also occur in more sophisticated preconditioning tech- 
niques for coupled systems of PDEs (see [33] for further details). The unknown 
function p in (4.1 ) is the pressure; a £ [L°°(fl)] 3x3 is the diffusion tensor and it usu- 



ally depends on the permeability and viscosity, etc.; c £ L°°(r2) is nonnegative and 
proportional to the inverse of the time step size in an implicit temporal scheme (such 
as the Backward Euler method). We assume that dtt has two non-overlapping parts, 



and Tff, such that To U Tn — dfl. The model problem (4.1 ) is completed by the 
boundary conditions p = go on and (aVp) n — g^ on T^r, where n is the outward 
unit normal vector of 90, and where go and g^ are given. 



The difficulties that arises from (4.1 ) in regard to solving the linear system (2.1 ) 
are mainly due to the complex geometry and complicated physical properties of the 
reservoir, which often result in a heterogeneous diffusion tensor with large jumps, and 
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in distorted, degenerated, and/or non-matching meshes with faults and pinch-outs. 
Figure |4.1| shows an example of the computational domain and the permeability of 
this reservoir in one horizontal layer. 



We use the mixed-hybrid finite element method to discretize (4.1). Due to its 
local conservation property and intrinsic and accurate approximation of the flux u := 
— aVp, the mixed-hybrid method is preferred in reservoir simulations. Let Qh — {T} 
be a triangulation of fi, and let f\ be the set of boundaries of T in f\ with the 
decomposition 

it^reivredfi}, r£ = r fc \r£ 

We define the finite dimensional spaces as 

V h := {v G [L 2 (n)] d : v\ T G V h (T), VT € il h } 

Qh := {q G L 2 (n) : q\ T G P (T), VT G n h } 

A h := {p G L 2 {T h ) : M | r G P (T), Vr G T°; M | r = 0, VF e T 9 h }, 

where Vh(T) is the Kuznetsov-Repin element on the polyhedral elements [IT]- This 
element is based on the Raviart-Thomas element on the local conforming tetrahedral 
partitioning of T. Pk(T) denotes the set of polynomials of a degree not greater than k, 
k > 0. Vi l: Qh, and A h are used to approximate the flux u; the pressure p associated 
with elements in f^; and the Lagrange multipliers Aft are associated with faces in fl^. 



The mixed-hybrid finite element formulation for (4.1) can be written as follows: 
Find (u h ,p h , Aft) G V h x Q h x K h , u h n = -gN.h on T h , where g N , h is an appropriate 
approximation of g^, such that 

{a^Uh, v) -^2{(p h , div v) T - {\ h ,v -n T ) dT } = - g D (v-n)ds, VveV h 
rp J r D 

- ^(div u h , q)r ~ (cph, q) = - J fqdx, Vg G Q h 

^2{p,u h ■ n T )dT = 0, V/i G A,,. 



T 

where denotes the outward unit normal to dT. 
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5. Numerical Experiments. It is easy to see that the resulting discrete linear 
system of the mixed-hybrid finite element method in the previous section has the 
following matrix form: 








3- 






\o 



Notice that the matrix D is block diagonal with each block corresponding to the flux 



unknowns in one element. Hence, we can easily invert D and obtain (2.1) with 



(5.1) A 



BD^B -4 




and / 



BD" 1 /! - h 
. CD- 1 /! , 



where it is well-known that A is SPD. 

In the numerical experiments, performed on a Dell Precision desktop computer, 
we solve the above linear system by PCG with different preconditioners. We pick 



eight problems from the real petroleum reservoir data. Table 5.1 gives the degrees of 
freedom (DOFs) and the number of non-zeros (NNZ) for the test problems. Here, the 
difference between Model 1 (3) and 2 (4) is that the permeabilities in Model 1 (3) are 



homogeneous whereas the permeabilities in Model 2 (4) have large jumps. Figure 4.1 
shows the computational domain of Models 3 and 4, and the highly heterogenous 
permeabilities used in Model 4. Models 5-8 are test problems that are relatively large 
in size. The computational domain of Models 5 and 6 are the same, but the physical 
properties, such as porosity and permeability, are different. Models 7 and 8 share the 
same computational domain, but the permeabilities used in Model 7 are from real 
reservoir data and the permeabilities in Model 8 are artificially adjusted in order to 
make the resulting linear system very difficult to solve. 

Table 5.1 

Degree of freedom (DOF) and number of non-zeros (NNZ) of the model problems. 



Test Problem Set 1 


Model 1 


Model 2 


Model 3 


Model 4 


DOF 


156,036 


156,036 


287,553 


287,553 


NNZ 


1,620,356 


1,620,356 


3,055,173 


3,055,173 


Test Problem Set 2 


Model 5 


Model 6 


Model 7 


Model 8 


DOF 


1,291,672 


1,291,672 


4,047,283 


4,047,283 


NNZ 


13,930,202 


13,930,202 


45,067,789 


45,067,789 



Table 5.2 compares different preconditioners for the model problems. The CPU 
time listed in Table |5.2| is the total CPU time which includes both the setup and 
solver times. We consider only the ILU(0) preconditioner due to memory-usage con- 
siderations. As we use highly unstructured meshes, the amount of fill-in for obtaining 
ILU(fc) factorizations is not predictable for k > 0. Here, the AMG(GS) refers to the 
V-cycle MG method with standard GS smoother on each level. The AMG(ILU(0)) 
refers to the V-cycle MG method with an ILU(0) smoother on the finest level and a 
GS smoother on the other levels. In both cases, two cycles are applied in order to 
define the preconditioner. B co is defined in (2.2 1 with one V-cycle MG, a GS smoother 



on each level serves as the smoother S, and ILU(0) serves as the preconditioner B. 
B is its additive version defined in (2.8) with the special choices of J = 1, W\ = V, 
III = I and A^ 1 is substituted by the preconditioner B. 
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Table 5.2 

Comparison of preconditioners for model problems (stopping criteria: relative residual less than 
10~ 10 ). In the table, '—' means the method does not converge within 10000 iterations, #Iter is the 
number of iterations, and the unit for CPU time is second. 



lest rroblcm oet l 


Model 1 


Model 2 


Model 3 


Model 4 




#Itcr 


CPU 


#Itcr 


CPU 


#Iter 


CPU 


#Itcr 


CPU 




234 


2.85 


3034 


32.91 


291 


6.56 


4939 


102.13 


AMG(GS) 


135 


13.16 


136 


13.48 


154 


31.70 


138 


25.78 


AMG(ILU(0)) 


5364 


631.72 


1058 


115.71 






4849 


1024.21 


B (Additive) 


44 


6.18 


43 


5.53 


44 


11.40 


49 


11.33 


B co (Multiplicative) 


23 


4.04 


25 


3.92 


22 


7.36 


28 


7.78 


Test Problem Set 2 


Model 5 


Model 6 


Model 7 


Model 8 


Preconditioner 


#Itcr 


CPU 


#Itcr 


CPU 


#Itcr 


CPU 


#Itcr 


CPU 


ILU(O) 










838 


271.46 


2039 


646.77 


AMG(GS) 


30 


36.70 


30 


37.07 


19 


96.03 


22 


112.79 


AMG(ILU(0)) 


15 


32.16 


15 


31.58 


11 


96.72 


11 


103.52 


B (Additive) 


32 


39.72 


32 


39.87 


26 


100.78 


25 


104.29 


B co (Multiplicative) 


19 


25.33 


19 


25.54 


14 


76.92 


14 


77.88 



The ILU(0) preconditioner works well for Models 1 and 3 because these two 
problems are relatively small in size and their permeabilities are homogeneous. We 
can see that the ILU(0) preconditioner does not work well for problems that are 
highly heterogeneous (Models 2 and 4). It is also not efficient for large-size problems 
(Models 5-8). Its performance is not predictable, and sometimes it may even break 
down (Models 5 and 6). AMG with a standard GS smoother works better than ILU(0) 
does, but still requires more than 100 iterations for Models 1-4. If we use ILU(0) as 
the smoother at the finest level, the resulting AMG preconditioner may not be SPD 
because ILU(0) might not converge. From the numerical results, we can see that for 
Models 1-4, using AMG with the ILU(0) smoother does not work as fast as using 
AMG with the GS smoother does. The former method may break down due to the 
lack of an SPD property (Model 3). However, if we combine AMG and ILU(0) as in 



Algorithm 2.1 the new defined preconditioner B co with AMG and ILU(0) gives the 
best performance in terms of CPU time; and, it performs efficiently and robustly with 
respect to the problem size and heterogeneity. The additive algorithm B also works 
efficiently and robustly. The numerical tests confirm our theoretical results and show 
the efficiency of our new approach for practical reservoir simulation problems. 



6. Concluding Remarks. In this paper, we discussed practical and efficient 
solvers for large sparse linear systems and their applications in petroleum reservoir 
simulation. We propose a simple and efficient framework for obtaining a new pre- 
conditioner by combining a convergent iterative method (such as an AMG method 
with a standard smoother) with an existing preconditioner (such as an ILU method). 
We used this approach to solve large linear systems arising from petroleum reser- 
voir simulation (with highly heterogeneous media and 3D unstructured grids). And, 
our numerical results show that the new preconditioners significantly improve the 
efficiency of the linear solvers and are robust with respect to the size and heterogene- 
ity of practical problems. This method has been implemented and employed in real 
reservoir simulation by ExxonMobil Upstream Company and RIPED, PetroChina. 



14 



Acknowledgement. The authors wish to thank Dr. Robert Falgout for many 
helpful discussions. 

REFERENCES 

[1] M. Benzi. Preconditioning techniques for large linear systems: A survey. Journal of Computa- 
tional Physics, 182(2) :418-477, 2002. 
[2] J. Bramble. Multigrid methods. Chapman & Hall/CRC, 1993. 

[3] A. Brandt. General highly accurate algebraic coarsening. Electronic Transactions on Numerical 

Analysis, 10(l):l-20, 2000. 
[4] A. Brandt, S. McCormick, and J. Ruge. Algebraic multigrid (AMG) for sparse matrix equations. 

In Sparsity and its applications (Loughborough, 1983), pages 257-284. Cambridge Univ. 

Press, Cambridge, 1985. 

[5] M. Brezina, R. Falgout, S. MacLachlan, T. Manteuffel, S. McCormick, and J. Ruge. Adaptive 
smoothed aggregation (aSA) multigrid. SIAM Rev., 47(2):317-346, 2005. 

[6] N. I. Bulccv. A numerical method for the solution of two-dimensional and three-dimensional 
equations of diffusion. Math. Sb., 51:227-238, 1960. 

[7] T. F. Chan and H. A. Van der Vorst. Approximate and incomplete factorizations. In Parallel 
numerical algorithms (Hampton, VA, 1994), volume 4 of ICASE/LaRC Interdiscip. Ser. 
Sci. Eng., pages 167-202. Kluwcr Acad. Publ., Dordrecht, 1997. 

[8] A. J. Cleary, R. D. Falgout, V. E. Henson, J. E. Jones, T. A. Manteuffel, S. F. McCormick, 
G. N. Miranda, and J. W. Ruge. Robustness and scalability of algebraic multigrid. SIAM 
J. Sci. Comput., 21(5):1886-1908 (electronic), 2000. Iterative methods for solving systems 
of algebraic equations (Copper Mountain, CO, 1998). 

[9] T. Clees and L. Ganzcr. An efficient algebraic multigrid solver strategy for adaptive implicit 
methods in oil-reservoir simulation. SPE Journal, 15(3):670-681, 2010. 
[10] O. Dubois, I. Mishev, and L. Zikatanov. Energy minimizing bases for efficient multiscale mod- 
cling and linear solvers in reservoir simulation. In SPE Reservoir Simulation Symposium, 
SPE 119137, 2009. 

[11] D. J. Evans. The use of pre-conditioning in iterative methods for solving linear equations with 
symmetric positive definite matrices. J. Inst. Math. Appl., 4:295-314, 1968. 

[12] R. Falgout. An introduction to algebraic multigrid. Computing in Science and Engineering, 
8(6):24, 2006. 

[13] W. Hackbusch. Multi-grid methods and applications, volume 4. Springer Verlag, 1985. 

[14] X. Hu, W. Liu, G. Qin, J. Xu, Y. Yan, and C. Zhang. Development of a fast auxiliary subspacc 
pre-conditioner for numerical reservoir simulators. In SPE Reservoir Characterization and 
Simulation Conference, SPE 148388, 2011. 

[15] Y. Jiang and H. Tchclepi. Scalable multi-stage linear solver for coupled systems of multi- 
segment wells and complex reservoir models. In SPE Reservoir Simulation Symposium. 
SPE 119175, 2009. 

[16] J. E. Jones and P. S. Vassilevski. AMGe based on element agglomeration. SIAM J. Sci. 

Comput., 23(1):109-133 (electronic), 2001. 
[17] Y. Kuznetsov and S. Repin. New mixed finite element method on polygonal and polyhedral 

meshes. Russian J. Numer. Anal. Math. Modelling, 18(3):261-278, 2003. 
[18] I. Lashuk and P. S. Vassilevski. On some versions of the element agglomeration AMGe method. 

Numer. Linear Algebra Appl, 15(7):595-620, 2008. 
[19] S. V. Nepomnyaschikh. Decomposition and fictitious domains methods for elliptic boundary 

value problems. In Fifth International Symposium on Domain Decomposition Methods for 

Partial Differential Equations (Norfolk, VA, 1991), pages 62-72. SIAM, Philadelphia, PA, 

1992. 

[20] J. W. Ruge and K. Stuben. Algebraic multigrid. In Multigrid methods, volume 3 of Frontiers 

Appl. Math., pages 73-130. SIAM, Philadelphia, PA, 1987. 
[21] Y. Saad. Iterative methods for sparse linear systems. Society for Industrial and Applied 

Mathematics, Philadelphia, PA, second edition, 2003. 
[22] K. Stuben. Algebraic multigrid (AMG): experiences and comparisons. Appl. Math. Comput., 

13(3-4):419-451, 1983. 

[23] K. Stuben, P. Delaney, and S. Chmakov. Algebraic multigrid (AMG) for ground water flow and 
oil reservoir simulation. In Proceedings of the Conference MODFLOW and More, pages 
17-19, 2003. 

[24] U. Trottenberg and T. Clees. Multigrid Software for Industrial Applications-From MG00 to 
SAMG. Numerical Flow Simulation I: Cnrs-Dfg Collaborative Research Programme, Re- 



15 



suits 1996-1998, page 423, 2001. 
U. Trottenberg, C. Oosterlee, and A. Schullcr. Multigrid. Academic Press, 2001. 
P. Vanek, J. Mandel, and M. Brczina. Algebraic multigrid by smoothed aggregation for second 

and fourth order elliptic problems. Computing, 56(3):179— 196, 1996. International GAMM- 

Workshop on Multi-level Methods (Mcisdorf, 1994). 
R. S. Varga. Matrix iterative analysis. Prentice-Hall Inc., Englewood Cliffs, N.J., 1962. 
P. S. Vassilevski. Multilevel block factorization preconditioners. Springer, New York, 2008. 
C. Wagner. Introduction to algebraic multigrid. Course Notes of an Algebraic Multigrid Course, 

University of Heidelberg, 1999. 
J. Watts III. A conjugate gradient-truncated direct method for the iterative solution of the 

reservoir simulation pressure equation. Old SPE Journal, 21(3):345-353, 1981. 
J. Xu. Iterative methods by space decomposition and subspacc correction. SIAM Rev., 

34(4):581-613, 1992. 

J. Xu. The auxiliary space method and optimal multigrid preconditioning techniques for un- 
structured grids. Computing, 56(3):215-235, 1996. 

J. Xu. Fast Poisson-based solvers for linear and nonlinear PDEs. In R. Bhatia, editor, Pro- 
ceedings of the International Congress of Mathematics, volume 4, pages 2886-2912. World 
Scientific, 2010. 



